Model for reversible colloidal gelation 
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We report a numerical study, covering a wide range of packing fraction and temperature T, for a 
system of particles interacting via a square well potential supplemented by an additional constraint 
on the maximum number nmax of bonded interactions. We show that when rimax < 6, the liquid-gas 
coexistence region shrinks, giving access to regions of low where dynamics can be followed down 
to low T without an intervening phase separation. We characterize these arrested states at low 
densities (gel states) in terms of structure and dynamical slowing down, pointing out features which 
are very different from the standard glassy states observed at high tf) values. 

PACS numbers; 82.70.Gg, 82.70.Dd, 61.20.Lc 



Extensive investigations have recently focused on slow 
dynamics in colloidal systems, triggered by the exper- 
imental discovery of multiple mechanisms leading to a 
disordered arrested state. At high packing fraction 0, ar- 
rest takes place via a glass transition process, which can 
be driven by jamming, as in hard-sphere (HS) systems, 
or by attractive bonding between colloidalparticles, gen- 
erating the so-called "attractive glass" 13 • At low 0, 
arrest takes place via particle clustering!^ and formation 
of an amorphous state of matter — named a gel — able 
to support weak stresses Clustering and gelation have 
no counterpart in atomic or molecular systems, being in- 
duced by the presence of attractive interactions between 
particles, with a range of interaction much smaller than 
the particle size. 

Recently, several studies have addressed the ques- 
tion of the routes to the gel state in colloidal sys- 
tems |lllllli|ll|ll|l2. In most models proposed so 
far, phase separation or microphase separation provides 
the initial step of the gelation process. For short-ranged 
attractive colloids, at low cf) and temperature T, the phase 
diagram is characterized by a flat phase-coexistence curve 
(e.g. see Fig. A quench inside the two-phase region in- 
duces, via spinodal decomposition, a separation into col- 
loid rich (liquid) and colloid poor (gas) phases. However, 
in appropriate conditions, dynamical arrest in the denser 
region, of the attractive glass type, intervenes by freezing 
the pattern generated during the coarsening process. The 
time to arrest and the structure of the formed gel depends 
on <j) and on the interaction strength. The connection be- 
tween gelation and phase separation 0, 0, 0, 0, 
is supported by the experimental observation of a peak 



in the scattered intensity at small wavevectors'4f, 18]. In 
the case of phase separation induced gelation, the slowing 
down of the dynamics associated with the formation of an 
arrested state can not be continuously followed through 
equilibrium states, since gelation takes place only after 
a quench inside the coexistence region. This is different 
from the slowing down of the dynamics on approaching 



the glass transition at high densities, for which the <f) and 
T dependence of the slow dynamics can be characterized 
in terms of (metastable) equilibrium states in a repro- 
ducible and reversible way. 

In principle, gelation could also occur in the absence 
of any phase separation, if the liquid-gas separation is 
suppressed or if the gel line is encountered before the 
liquid-gas coexistence locus'^, preventing phase separa- 
tion. In these cases, the approach to the gel state could be 
in principle followed in equilibrium and reversibly. Some 
experimental groupsfisj do indeed favor this interpreta- 
tion, explaining the increase in the scattered intensity as 
intrinsically due to the inhomogeneities associated with 
the gel structure, posing a challenge to theoreticians to 
develop a model where arrest at low (/> is observed in the 
absence of a phase separation mechanism and in which 
reversible physical gelation can be studied in equilibrium. 

In this Letter we introduce and numerically study a 
simple model showing that dynamical arrest at low (as 
low as 0.2) can indeed be generated in the absence of 
phase separation. We follow, via extensive simulations, 
the evolution, in equilibrium, of the density autocorre- 
lation functions to point out the differences in the dy- 
namics of a gel and of a glass, and to provide a way of 
discriminating between the two phenomena. 

We study a modification of the short-range square well 
(SW) potential, by adding a constraint which limits the 
maximum number of bonds rimax that can be formed 
by the particles Without the constraint, each par- 
ticle could in principle interact simultaneously with the 
maximum number of neighbours allowed by geometrical 
packing (rimax = 12). The constraint switches off the 
attractive well potential when any of the two interact- 
ing particles has already rimax bonded neighbors. In this 
way, the energy difference between particles located in 
the interior of an aggregate and particles at the aggre- 
gate surface can be decreased and even suppressed. As 
a result, the surface tension decreases poj. the driving 
force for phase separation is significantly reduced and 
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open structures are favored. We find that for rimax < 6, 
the system can access the region usually dominated 
by phase separation and experience a dynamical slowing 
down by several orders of magnitude, thus entering the 
gel regime. Remarkably, the system undergoes the fluid- 
gel transition in equilibrium, so that the process is fully 
reversible. 

We perform event-driven molecular dynamics (MD) 
simulations of particles interacting via a maximum va- 
lency model^l], i.e. a SW potential, where particles 
can form a maximum number of bonds rimax- We fix 
the depth uq = 1, and the width A/(cr -I- A) = 0.03, 
with (7 = 1 being the particle hard-core diameter. The 
square well form of the potential unambiguously defines 
bonded particles when particle centers lie within a dis- 
tance a < r < a + A from each other. When a parti- 
cle is already bonded to rtmax neighbors, the well inter- 
action with other particles is switched off leaving only 
the HS potential. We simulate a large system containing 

= 10^ particles of mass m = 1 to minimize finite size 
effects. Temperature is measured in units of uq, time t in 
units of a{ra/uf)Y/'^ . For all simulated state points, we 
first equilibrate the system at constant T until the po- 
tential energy and pressure P of the system have reached 
a steady state, and particles have diffused several cr's on 
average. Then, a production run is performed at constant 
energy. An average over typically 100 different realiza- 
tions is done to gather statistics. We focus on the cases 
Hmax = 3, 4 and 5, since nmax = 6 behaves similarly to 
the unconstrained case as far as liquid-gas phase separa- 
tion is concerned. We perform MD simulations to effi- 
ciently propagate the system in time, even if Brownian 
dynamics (BD) would give a more realistic description 
of the short time dynamics. As far as slow dynamics is 
concerned, MD and BD are equivalent |23.l23|. 

We report in Fig. ^ the evolution of the spinodal line 
for different nmax in the (0, r)-plane. We estimate the 
spinodal line by bracketing it with the last stable state 
point and the first phase separating state point along 
each isochore. The last stable state is characterized by a 
value of the structure factor at low q of the order of 10. 
We confirm the location of the spinodal by detecting the 
vanishing of the derivative of P{V) along isotherms. The 
unstable area in the (i^, T)-plane shrinks on decreasing 
^max, showing that the additional constraint opens up a 
significant portion of phase space, where the system can 
be studied in equilibrium one-phase conditions. The (p 
at which phase separation is not present, at all T, can 
be as low as « 0.20 for rtmax — 3 0|, ~ 0.30 for 
^max — 4 and (j) « 0.35 for rimax = 5. Fig. ^ also shows 
the percolation lines, calculated as the loci where 50% 
of independent configurations are characterized by the 
presence of a spanning cluster of bonded particles. The 
percolation loci also shift to lower T on decreasing nmax, 
always ending in the spinodal on the low (j) side. 

As shown by the data in Fig^ the addition of the con- 
straint on the maximum number of bonds, by suppressing 
the phase separation, makes it possible to study the dy- 




FIG. 1: Spinodal (dashed) and percolation (full) loci (with 
reported error bars) for Wmax = 3 (circles), 4 (squares), 5 
(triangles) and 12 (no symbols). Lines are guide to the eye. 



namics of the model at very low T, where the lifetime of 
the interparticle bond increases, stabilizing for longer and 
longer time intervals the percolating network [25j. When 
the bond lifetime becomes of the same order as the obser- 
vation time, the system will behave as a disordered solid. 
It is worth stressing that in the present model there is no 
thermodynamic transition associated with the onset of a 
gel phase p^. 

To quantify these propositions we study in the follow- 
ing the dynamics for all T in the </> region where the 
system is in a single phase. 
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FIG. 2: Arrhenius plot of the bond lifetime tb for different 
values of rimax and (p. The vertical line indicates the lowest 
accessible T in the one-phase region for the unconstrained 
(jimax = 12) case. Dashed lines are Arrhenius fits. 

Fig. 121 shows the T-dependence of the bond lifetime 
Tb, defined as the time at which the bond autocorrela- 
tion function decays to 0.1, along different isochores for 
"max = 3, 4 and 5. The bond lifetime follows an Arrhe- 
nius behavior at low T, with an activation energy slightly 
increasing with nmax, of the order of the bond energy, 
suggesting that breaking of one single bond is the ele- 
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FIG. 3: (r^> for n,nax = 3 and = 0.20 

as a function of T. From left to right, Ts are: 
1.0,0.5,0.3,0.2,0.15,0.125,0.1,0.09,0.08. The arrow signals 
the typical plateau value for the HS glass. 



mentary process. We do not observe any significant (j) 
dependence, suggesting that the bond- breaking depends 
mostly on T|28j as well as on nmax- AH curves, for all re- 
ported n„iax and 0, superimpose onto each other for high 
T, recovering the HS limit. The vertical line in Fig. |21 
indicates the lowest T which can be reached before en- 
countering the coexistence region for the unconstrained 
case. It is interesting to note that the introduction of a 
small rimax makes it possible to explore states with dy- 
namics at least 4 orders of magnitude slower than with- 
out the constraint, allowing for an approach to arrested 
states from equilibrium conditions. 

Fig. O shows the T dependence of the mean squared 
displacement (r^). A T-independent plateau develops on 
cooling, indicating the presence of a localization length of 
order cr, much larger than the corresponding value typi- 
cal of glass forming systems ( (r^) ~ O.Olcr^). The local- 
ization length decreases on increasing (jj. Caging in the 
gel is thus induced neither by the bond length (A) nor 
by the nearest neighbor distance and hence it is signifi- 
cantly different from the case of attractive and repulsive 
glasses 0, 13 • This feature signals the crucial role of con- 
nectivity in the gel transition. 

To quantify the dynamical slowing down on approach- 
ing the gel transition, we study the normalized density 
autocorrelation functions = {pq{t)p-q{0)) /S{q). 

Fig0| shows the behavior of (pqit) on decreasing T. A 
striking dependence on q is observed, a feature missing 
in the slow dynamics close to the glass transition. If we 
focus, as typically done in glass transition studies, on the 
q value of the first peak of S{q), corresponding to the 
inverse average nearest-neighbor distance, we observe no 
sign of a plateau in (pqit), within the precision of our 
calculations (left panel of Fig. However, at smaller 
q values, a clear plateau, named fq, in analogy with its 
counterpart for the glass transition, where it is usually 
named non-ergodicity factor, emerges at a T that varies 
with q. This behavior is profoundly different from what 
is observed in a standard glass transition, but closely 
resembles what is observed at a percolation transition in 
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FIG. 4: Density autocorrelation functions (j>qi't) for f^max ~ 3 
and (fi — 0.20 as a function of T (sequence of values as in 
FigEJi respectively for qa ~ 2n (left) and gcr ~ 1 (right). 
Note the difference in T dependence of in the two cases. 
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FIG. 5: Left panel: Plateau height fq for rimax = 3 and (j> = 
0.20 as a function of T, for different values of qa. At a given 
T, fq at small q. At larger q, though qa is still smaller 
than 27r (nearest neighbor length scale), fq ~ within the 
accuracy of the data. Right panel: at T = 0.09, (f> = 0.20, 
^imax = 3, compared with fq of attractive and IfS glass ||30||. 



the presence of chemical (infinite lifetime) bonds |12ll29| . 
We find that the T-dependence of the relaxation time Ta 
(defined as the time at which 4>q{t) reaches the value 0.1) 
crosses over from ^ 1/VT at nearest neighbour dis- 
tances to Ta ^ e^l^ at large length scales. For small 9, 
Ta is coupled to the bond lifetime, while at larger 5, 
is controlled by the microscopic time scale (proportional 
to the thermal velocity). Fig. shows the T-dependence 
of the plateau value fq (calculated as the amplitude of a 
stretched exponential fit to the slow relaxation) for vari- 
ous g's (left panel), and compares its q dep endence with 
/q's typical of HS and attractive glasses |^ (right panel). 
While in the case of glasses (either HS or attractive) fq is 
significantly different from zero at all physically relevant 
q values, in the present case fq is significantly different 
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from zero only at very small wavevectors. It is also in- 
teresting to observe that, while in the case of glasses fq 
jumps from zero to a finite value at the glass transition, 
in the present case fq appears to grow smoothly from 
zero. 

In summary, we have proposed an off-lattice model 
which allows the study of dynamical arrest at low 
(physical gelation) in the absence of macro or micro phase 
separation, effectively decoupling the effects of phase sep- 
aration from the dynamical slowing down. The process 
of arrest can be followed in equilibrium. This condition 
offers the possibility of studying the behavior of the den- 
sity correlators and of the mean squared displacement 



close to the gel transition, which is strongly coupled to 
the bond dynamics. We discover significant differences 
between glasses and gels in terms of the decay of the den- 
sity correlation functions and of the localization length, 
which could be studied experimentally to discriminate 
between the phase separation route and the single phase 
reversible gelation case. 
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